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SOLUTION  OF  PARTIAL  DIFFERENTIAL  EQUATIONS 
BY  MEANS  OF  FINITE  ELEMENTS 

AN  INTRODUCTORY  SKETCH 
J.  Barkley  Rosser 

1.  The  mathematical  framework.  Engineers  have  long  found  it 
advisable  to  analyse  the  stresses  and  strains  in  structures  that  they 
design  or  construct.  Not  uncommonly,  these  structures  are  composed 
of  many  parts.  By  appealing  to  such  principles  as  virtual  work,  equations 
can  be  written  relating  the  stresses  and  strains  of  one  part  to  those  of 
adjacent  parts.  Trom  this  there  results  a set  of  simultaneous  linear 
equations,  which  can  be  solved  to  get  a complete  picture  of  all  the 
stresses  and  strains.  For  a structure  with  very  many  parts,  the  number 
of  simultaneous  equations  will  be  very  great.  However,  since  most 
parts  are  adjacent  to  only  a few  other  parts,  most  coefficients  are  zero. 
Hence,  by  ordering  the  equations  judiciously,  solutions  can  be  obtained 
fairly  quickly.  Because  the  parts  are  separate  and  finite  elements,  the 
procedure  became  known  as  the  method  of  finite  elements.  A short  list  of 
references  to  books  on  this  subject  is  given  at  the  end  of  the  report.  The 
most  comprehensive  one  is  that  by  Zienkiewicz,  but  others  may  be  more 
suitable  for  inexperienced  readers. 
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The  method  had  such  success  that  it  was  adapted  to  many 


situations  in  which  there  were  finite  elements  only  by  arbitrary  definition. 

One  such  area  is  the  approximate  numerical  solution  of  differential  equations. 
We  shall  give  an  introductory  sketch  of  the  method  as  applied  to  certain 


partial  differential  equations.  Specifically,  consider  the  equation 


a au 
ax  p ax 


a , au 

— (q  — ) + ru 

ay  ay 


where  u,  f,  p,  q,  and  r are  functions  of  x and  y.  We  interpret  L 
as  that  operator  on  u,  depending  on  p,  q,  and  r,  such  that  Lu  is  the 
left  side  of  (1.1).  The  specific  problem  is  to  find  a u which  satisfies 
(1.1)  inside  a given  region  £2  and  satisfies  some  given  conditions  on 
the  boundary  r (£2)  of  £2.  For  example,  we  could  ask  that  u satisfy 
(1.11  inside  the  square  0 < x < 1,  0 < y < 1,  and  ask  further  that  u 0 
on  the  sides  of  the  square. 

We  are  not  here  seeking  great  generality.  We  shall  consider 
regions  £2  which  are  fairly  simple,  and  have  well  behaved  boundaries. 
Also,  we  shall  ask  that  each  of 

2 2 
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be  continuous.  This  certainly  assures  that  the  left  side  of  (1.1)  makes 
sense;  indeed  is  continuous. 
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Wa  define  the  inner  product  (u,  v)  by 


(1.  Z) 


(u,  v)  J J uv  dx  dy  . 


12 


Theorem  1.1.  The  equation  (1.1)  holds  in  12  if  and  only  if 


(1.  3) 


(Lu-f,  v)2  = 0 


for  all  v which  are  continuous  on  12  and  equal  to  0 on  T (12). 

Proof.  if  (1.1)  holds,  then  obviously  (1.3)  holds.  Now  assume 


that  (1.  3)  holds  as  stated.  If  (1.  3)  should  hold  for  v = Lu-f,  giving 


(1.4) 


J J (Lu-f)Zdx  dy  = 0, 


then  obviously  Lu  - f would  have  to  be  zero  all  over  12;  that  is,  (1.1) 
holds.  However,  Lu-f  is  likely  not  zero  on  r(i2).  What  we  do  is 
take  v Lu  - f except  very  close  to  the  boundary  of  12,  and  then  bend 
it  to  be  0 at  the  boundary.  Then  the  left  sides  of  (1.  3)  and  (1.4)  will  be 
nearly  equal.  By  restricting  the  region  where  we  bend  Lu-f  to  be 
closer  and  closer  to  the  boundary,  we  make  the  left  sides  of  (1. 3)  and 
(1.4)  as  nearly  equal  as  we  like.  So  (1.  4)  must  hold. 

We  define  the  operator  M(u,  v)  by 


U.  51  M(u,  »)  = //  (P  g fjf  + <1  JZ  S + rUV)dX  dy- 


12 


?y  By 
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1 


Theorem  1 . 1.  If  v 0 on  r'  ) and  3v/9x  and  3v/3y  are 
continuous,  then 


(1.6) 


M(u,  v)  = (Lu,  v)  . 


Proof.  In  texts  on  advanced  calculus,  Green's  theorem  is  stated 
in  such  forms  as 


(1.7) 


Take 


/ (Adx  + Bdyt  = J J (“  - — )dx  dy 


1(S2) 


Si 


9x  ?y 


A = -vq 


ay 


D ?'U 

B = vp  — 

f>x 


Then  the  left  side  of  (1.7)  is  zero,  because  v = 0 on  r(^).  So  we  get 

J J { P + q f ~ Idx  dy  = - f f v{-^-  (p  — ) + -^-(q  — )}dx  dy. 

J J 3x  3x  3y  3y  »v  av  av  m av  1 


3x  ax  3y  ?y 


Si  Si 

From  this,  (1.6)  easily  follows. 

We  say  that  u is  a generalized  solution  of  our  problem  if  it  has  continuous 
first  derivatives  and  satisfies  the  given  conditions  on  r (Q)  and  if 


(1.8) 


Mfu,  v)  = (f,  v)2 


for  every  v which  is  0 on  r (n)  and  has  continuous  first  derivatives. 
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If  u is  a generalized  solution,  then  by  Theorem  1.2,  we  have 
(Lu  - f,  v)  0 for  every  v which  is  0 on  r (ft)  and  has  continuous  first 

L. 

derivatives.  However,  we  cannot  use  Theorem  1.1  directly  to  conclude  that 
u satisfies  (1.1),  because  Theorem  1.1  requires  that  (1.  3)  hold  for  still  more 
v's.  However,  unless  we  are  dealing  with  some  very  poorly  behaved  functions, 
it  will  usually  be  the  case  that  if  u is  a generalized  solution,  it  will  indeed 
satisfy  (1.1). 

By  reasoning  as  in-th®  proof  of  Theorem  1.1,  we  conclude  that  if  u is 
a generalized  solution,  then  (1.8)  holds  for  every  v with  continuous  first 

derivatives . 

Theorem  1.3.  If  p,  q,  and  r are  non-negative,  then  if  u 
is  a generalized  solution,  it  is  a w with  continuous 
first  derivatives  taking  the  given  values  on  F (ft)  that  minimizes 

(1.9)  M(w,  w)  - 2(f , w)  . 

Proof.  Let  w minimize  (1.9).  Take  v with  continuous  first 
derivatives  and  equal  to  £ on  r(ft).  Then  for  real  a , w + av  will  have 
continuous  first  derivatives  and  take  the  given  values  on  F(ft),  So 

M(w  + av,  w + av)  - 2(f,  w + av) 

must  have  a minimum  at  a = 0.  So  its  derivative  must  be  zero  at  a - 0. 


Expanding  it  gives 


M(w,  wi  + 2uM(w,  v)  + aZMlv,  v)  - 2(f,  w)_,  - 2o(f,  . 

Differentiating  with  respect  to  a and  setting  0 gives 

Miw,  v)  - (f,  v)  . 

As  this  holds  for  all  suitable  v,  we  conclude  that  w is  a generalized 
solution.  Alternatively,  let  u be  a generalized  solution,  and  let  v be 
another  function  with  continuous  first  derivatives.  Consider 

Q {M(v,  v)  - 2 ( f , v)  1 - (M(u,  u)  - 2(f , u^l  . 

Taking  v ~ u in  (1.8)  gives 

Q - Mi'v,  v)  + M(u,  u)  - 2(f,  v)  . 

Using  (1.8)  in  this  gives 

Q = M(v,  v)  + M(u,  u)  - 2M(u,  v) 


= M(u-v,  u-v) 


-SJ  fp< 


a(u-v), 

ax 


2 


,a(u-v)  2 , , ,2  , , 

+ q( -)  + r(u-v)  Idxdy. 

ay 


As  p,  q,  and  r are  non-negative,  we  get  Q > 0.  So  u does  produce 
minimum. 

With  slightly  stronger  conditions  on  p,  q,  and  r,  we  easily 
conclude  that  u is  the  unique  w that  minimizes  (1.9). 
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This  leads  to  the  Rayleigh -Ritz  procedure  for  approximating  a 
generalized  solution.  Let  u^,  u^,  u^  all  have  continuous  first 

derivatives.  Let  uQ  take  the  given  values  on  r($2).  Let  u^,  u2,  ... 
all  be  zero  on  I'(S2).  Let  c^  c^,  cn  be  real  parameters . Then 


(1.10) 


n 

u0  + L CkUk 

k 1 K K 


has  continuous  first  derivatives  and  takes  the  given  values  on  I' (12). 

W'i  substitute  (1.10)  for  w in  (1.9),  and  choose  the  c's  so  as  to  minimize 
the  resulting  expression.  To  make  this  a minimum,  the  partial  derivatives 
with  respect  to  the  various  c^'s  must  all  be  zero.  Taking  the  partial  with 
respect  to  gives 


n au.  au  ?u  ?u 

l \ JJ  <»  -d  iz 

k 1 12 


+ q — — — + ru.u,  Idx  dy 
3x  dx  py  Py  ) k 


(1.11) 


pu  au.  au  au. 

JJ  (P—  T^+q—  T^-  + runuJ  )dx  dy 


12 


ax  ax 


ay  ?y 


o j 


(f,  u. 


j 2 

This  is  a set  of  simultaneous  linear  equations  to  be  solved  for  the 
c's.  What  this  does  is  to  pick  out  the  best  approximation  (in  some  sense) 
of  the  form  (1. 10)  for  the  u in  question.  If  we  have  chosen  the  uk‘s 
cleverly,  so  that  there  is  a good  approximation  to  u of  the  form  (1.10), 
we  will  have  found  it. 
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A 


The  Galerkin  approach  proceeds  according  to  a different  principle. 

By  Theorem  1.1,  we  seek  a u that  satisfies  (1.3)  for  all  continuous  v equal 

to  0 on  I'  (12).  If  we  can  find  a u that  satisfies  (1.  3)  for  a judiciously 

chosen  set  of  v. , it  cannot  be  too  different  from  a u that  satisfies  (1.  3) 

for  all  v.  So  try  a u of  the  form  (l.lOj  and  let  us  require  that  it  satisfy 

(1.3)  for  v successively  taken  equal  to  v. , v_,  ...,  v . For  v there 

i 2 n j ’ 

results 


(1.12) 


y c,  (p 

Lj  k ox  F 


k = l 


8u, 


ax 


„ 0U„ 

JL  , 0 

8x  P ax 


_ JL 

(q 

k 

— -) 

+ ru.  , 

v L 

” ay 

ay 

k’ 

J 2 

_ JL 

(q 

auo 

) + ru„ 

, v.) 

ay 

ay 

0 

’ j 

- ,f'  hh- 


Again  we  have  a set  of  simultaneous  linear  equations  for  the  c . 

k 

The  resulting  value  of  (1.10)  not  only  satisfies  (1.3)  for  each  v , but  for 

any  linear  combination  of  the  v^ . If  the  v^'s  are  chosen  so  that  arbitrary 

functions  can  be  approximated  well  by  linear  combinations  of  the  v 's  we 

j ’ 

can  hope  that  we  have  come  close  to  satisfying  (1.  3)  for  arbitrary  v's. 

One  can  try  a different  attack  by  the  Galerkin  approach.  The  function 

u is  a generalized  solution  if  it  satisfies  (1.8)  for  arbitrary  v.  Let  us 

try  to  satisfy  (1.8)  for  v taken  successively  equal  to  v, , v_ v . 

1 Z’  ’ n 


For  v.  there  results 
J 
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It  will  be  noted  that  if  we  take  the  v's  the  same  as  the  u's,  this 
reduces  to  (1.11),  which  was  obtained  by  the  Rayleigh -Ritz  approach. 
However,  this  derivation  of  (1.11)  is  probably  preferable  because  in  deriving 
(1.11)  by  the  Rayleigh -Ritz  approach  we  worked  from  (1.9),  in  which  w was 
prescribed  to  take  certain  values  on  I'  (12) . For  a u which  is  prescribed  to 
satisfy  mixed  conditions  on  F(i2)  the  derivation  from  (1.8)  by  the  Galerkin 
approach  is  satisfactory. 

An  advantage  of  (1.13)  over  (1.12)  is  that  for  (1.12)  one  would  have 
to  choose  the  u's  to  be  twice  differentiable,  whereas  for  (1.13)  singly 
differentiable  u's  will  suffice.  This  will  be  an  important  consideration 
later  on. 

2.  Finite  elements  to  the  rescue.  We  now  face  up  to  the  real 
problem,  namely  the  judicious  choice  of  the  u's  (and  perhaps  the  v's) 
for  the  Rayleigh -Ritz  or  Galerkin  approaches.  It  would  be  helpful  if 
we  could  choose  them  so  that  the  integrations  inherent  in  (1.11),  (1.12), 
and  (1.13)  could  be  easily  performed.  This  is  not  absolutely  necessary, 
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since  we  could  grind  out  approximations  on  a computer  by  established 
quadrature  formulas.  It  would  be  even  more  helpful  if  we  could  arrange 
that  a large  fraction  of  the  coefficients  of  the  c's  would  be  zero.  This 
would  make  the  equations  much  easier  to  solve.  It  will  turn  out  that 
this  can  be  done. 

To  fix  our  ideas,  let  11  be  the  square  0 < x < 1,  0 < y < 1;  see 

Figure  1.  We  subdivide  this  into  triangles,  whose  vertices  form  a grid 


on  U.  With  each  grid  point  we  associate  an  area  of  support  and  a u^. 

The  area  of  support  consists  of  all  the  triangles  that  have  this  grid  point 

as  a vertex.  Then  for  each  grid  point  we  define  the  function  u^  as 

follows.  Outside  the  area  of  support,  is  identically  zero.  Inside 

the  area  of  support,  u^  is  a function  which  is  1 at  the  grid  point  and 

0 along  the  sides  of  the  triangles  opposite  the  grid  point.  One  very 

easy  way  to  do  this  is  to  have  consist  entirely  of  planes  joined 

together  along  the  sides  of  the  triangles.  Outside  the  area  of  support, 

u is  a horizontal  plane  at  zero  altitude.  For  each  triangle  in  the  area 
k 

of  support,  u is  the  plane  which  is  one  unit  high  at  the  grid  point  and 

fv 

zero  units  high  along  the  opposite  side  of  the  triangle.  That  is,  for  the 
area  of  support  consists  of  a pyramid  of  height  unity  at  the  grid  point, 
dropping  down  to  zero  at  the  edge  of  the  area  of  support,  while  u^  is 
identically  zero  outside  the  area  of  support. 

If  now  we  try  to  approximate  a function  u which  has  the  value 
at  the  k-th  grid  point, 

n 

(2.1)  ^ ckUk 

k=l 

is  obviously  an  approximating  function.  It  has  exactly  the  right  value 
at  each  grid  point.  Not  only  that,  but  if  the  function  u is  fairly  smooth, 
(2.1)  is  not  a bad  approximation  anywhere.  To  be  specific,  let  u have 
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continuous  derivatives  up  to  the  second  order,  and  let  be  the 

maximum  absolute  value  of  any  derivative,  first  or  second,  in  12. 
Then  the  difference  between  u and  (2.1)  is  at  most  4M  where 

h is  the  longest  side  of  any  triangle;  see  the  Lemma  on  p.  142  of 
Prenter. 

For  any  triangle  in  the  area  of  support,  it  is  easy  to  write 
the  equation  of  the  plane  that  constitutes  u^  above  that  triangle. 
Let  (xQ,  yQ),  (Xj,  y^),  and  (x^,  y^)  be  the  coordinates  of  the 
vertices  of  the  triangle.  We  seek  the  plane  which  is  one  unit  high 
at  (xR,  yQ)  and  zero  units  high  at  each  of  (x^,  y^)  and  (x^,  y^). 
Inside  the  triangle  we  would  have 


(xrx2)(yrY2)"(xrx2l(Y"Y2)+(x'xiHvrY2) 

,k  ° <'<1-x2l'V'i-y2l-l''rX2)(y0-y2H(Xo-x1)(yry2) 


The  denominator  is  twice  the  area  of  the  triangle,  and  so  is  not  zero. 
In  (1.8)  we  specified  that  v was  to  have  continuous  first 


derivatives.  But  no  u.  has  continuous  first  derivatives.  The  discon- 

k 

tinuities  are  not  bad,  just  along  the  sides  of  certain  triangles.  Else- 
where, each  u,  has  continuous  derivatives.  It  turns  out  that  this  is 
good  enough.  So  we  can  take  the  v.'s  in  (1.12)  or  (1.13)  to  be  the 
u^'s  as  defined  above.  Note  that  if  (1.12)  or  (1.13)  holds  for  all  the  v's, 
it  holds  for  any  linear  combination  of  them.  So  in  affect  we  are  taking 
v to  be  (2.1).  But,  as  we  saw,  any  v can  be  represented  fairly  closely 
by  (2.1).  So  we  are  close  to  assuring  that  (1.12)  or  (1.13)  holds  for  any  v. 
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It  seems  not  unreasonable  to  try  (2.1)  as  an  approximation  for  u. 

If  we  do  this,  use  of  (1.12)  will  not  be  practical,  since  it  requires  second 

derivatives  of  u,  whereas  for  the  u.  defined  above  even  the  first 

k 

derivatives  leave  something  to  be  desired.  However,  we  can  use  (1.13) 
perfectly  well.  So  we  use  (1.13),  taking  the  u,  's  and  v.'s  to  be  just 
the  u^'s  defined  above.  This  makes  (1.13)  the  same  as  (1.11),  except 
that  more  flexibility  is  permitted  for  the  behavior  of  u along  F (£2).  If 
the  values  of  u are  assigned  along  r (ft),  then  we  know  the  value  to  take 
for  Cj,  at  each  boundary  grid  point.  Thus  we  would  take  u^  to  be  the  sum  of  u^'s 
for  boundary  grid  points  multiplied  by  the  respective  boundary  values. 

In  view  of  (2.  2),  the  integrations  in  (1.13)  can  be  carried  out 
with  the  greatest  of  ease.  Also,  since  each  u^  is  identically  zero 
outside  the  area  of  support,  many  of  the  coefficients  of  the  c^  in  (1.13) 
will  be  zero;  if  we  have  a large  number  of  mesh  points  most  of  the 
coefficients  will  be  zero. 

If  we  knew  the  values  of  u at  each  grid  point,  we  could  take 
the  c,  to  be  these  values,  and  then  (2.1)  would  differ  by  less  than 

K 

2 

4M  h from  u at  any  point.  However,  the  c,  are  got  by  solving 

c.  K 

the  equations  (1.13).  The  resulting  c^  will  likely  not  be  exactly  the 
values  of  u at  the  k-th  grid  point.  Indeed,  there  might  be  some  question 
if  they  are  even  close.  For  the  special  case  p = 1,  q s 1,  and  r s 0, 
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Prenter  on  pp.  231- Z36  presents  a proof  derived  from  Friedrichs  that  the 


c are  close  to  the  values  of  u at  the  grid  points,  so  that  (2.1)  is 
k 

close  to  u.  In  fact  (see  the  top  of  p.  236  of  Prenter),  there  is  a constant 
M,  independent  of  how  the  grid  is  constructed,  so  that 

2 

(2.3)  Jj  (u-u  ) dx  dy  < ~ , 

S2  sin  0 

where  u stands  for  (2.1),  h is  the  longest  side  of  a triangle,  and 
is  the  smallest  angle  of  a triangle.  Hence,  by  going  to  a fine  enough  mesh, 
and  being  careful  not  to  use  triangles  with  small  angles,  one  can  contrive 
to  make  (2.1)  as  close  as  desired  to  the  solution.  Incidentally,  in  the 
third  displayed  formula  on  p.  233  of  Prenter,  one  should  put  K for 
""s/b-a  " . 

As  noted,  the  derivation  in  Prenter  assumes  p s 1,  q = 1,  and 
r = o.  These  requirements  can  be  relaxed  without  invalidating  the  result. 
However,  some  restraints  will  still  be  needed.  For  example,  in  Theorem 
1.3,  the  condition  that  p,  q,  and  r are  all  non-negative  was  invoked; 
if  one  wished  to  conclude  that  u is  the  UNIQUE  minimizing  function, 
slightly  more  was  required.  Of  course,  we  are  working  from  (1.13),  which 
did  not  depend  for  its  motivation  upon  Theorem  1.3.  However,  if  NO 
restraints  are  imposed  on  p,  q,  and  r,  there  may  fail  to  be  a solution 
for  (1.1).  In  any  case,  if  the  left  side  of  (1.1)  is  to  make  any  sense,  some 
sort  of  differentiability  seems  to  be  required  of  p and  q.  The  matter 
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seems  to  hinge  upon  whether  the  differential  operator  L is  "strongly 
coercive.  " The  short  discussion  on  p.  254  of  Prenter  seems  to  indicate 
that  if  p and  q have  continuous  first  partial  derivatives  and  are  positive 
and  bounded  over  12,  and  r is  continuous,  non-negative  and  bounded 
over  12,  then  one  may  expect  to  get  good  approximations  for  u by  solving 
(1.13).  Also,  the  errors  will  be  less  than  those  indicated  above  (see  p.  251 
of  Prenter).  If  p,  q,  and  r fail  to  satisfy  the  conditions  stated,  advice 
from  an  expert  should  be  sought.  Incidentally,  on  p.  254  of  Prenter,  read 
"Borman"  for  "Birman"  and  change  the  last  five  words  of  the  stated  theorem 
to  read:  "...  then  A is  strongly  coercive." 

It  seems  not  agreed  whether  in  this  approach  the  "finite  elements" 
are  the  triangles,  or  the  u^  based  on  the  grid  points,  or  something  else. 
Never  mind,  this  is  called  a method  of  finite  elements.  It  has  many 
variations.  For  one  thing,  it  can  be  used  for  solving  ordinary  differential 
equations;  see  pp.  201-227  of  Prenter. 

The  results  above  can  be  generalized  in  various  ways.  The  region 
12  may  be  subdivided  into  rectangles;  see  pp.  116-137  of  Prenter.  Or  one 
may  take  the  u^  to  consist  of  curved  surfaces.  A simple  way  to  do  this 
(discussed  on  pp.  144-153  of  Prenter)  proceeds  as  follows.  In  the  triangular 
grid  that  we  had  before,  introduce  new  grid  points  by  taking  one  in  each  side 
of  a triangle;  usually  the  midpoint  of  the  side  is  taken,  but  this  is  not 
obligatory.  As  before,  each  grid  point  is  to  have  an  area  of  support  and  a 
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uj.  • f or  vertices  of  triangles,  the  area  of  support  is  as  before.  Tor  each 

triangle  in  the  area  of  support,  u^_  is  given  as  the  product  of  two  formulas 

like  the  right  side  of  (2.  2)  (see  p.  147  of  Prenter)  and  so  is  a quadratic 

polynomial.  Tor  grid  points  on  the  sides  of  triangles,  the  area  of  support 

is  the  pair  of  adjacent  triangles  (only  one  triangle  if  the  grid  point  is  on 

the  boundary).  For  each  triangle  in  the  area  of  support,  u^  is  given  as 

the  product  of  two  formulas  like  the  right  side  of  (2.  2)  (see  p.  147  of 

Prenter).  Outside  the  area  of  support,  u^  is  identically  zero  in  all  cases. 

At  the  grid  point  itself,  u^  is  unity,  and  it  drops  off  to  zero  at  the  edge 

of  the  area  of  support;  indeed  this  u^  is  zero  at  all  other  grid  points. 

Then  (2.1)  again  serves  as  an  approximating  function.  If  u has  continuous 

derivatives  up  to  the  third  order,  and  M is  a bound  of  the  absolute  values 

3 

of  the  derivatives,  then  (see  p.  150  of  Prenter)  the  difference  between  u 
and  (2.1)  is  at  most 

8M3h3 
sin  0 ’ 

where  h is  the  longest  side  of  a triangle  and  0 is  the  least  angle 
of  a triangle. 

Although  these  u^'s  are  curved  surfaces,  and  make  (2.1)  fit 
u better  than  before,  they  still  have  discontinuous  derivatives  along 
sides  of  the  triangles.  Still  more  complicated  u^'s  have  been  devised, 
which  have  continuous  derivatives  everywhere.  This  increases  the 
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difficulty  of  calculation  quite  a bit,  and  it  is  not  at  all  clear  that  the 
gain  is  worth  the  extra  effort.  Many,  many  ingenious  ways  to  define  the 


u^'s  have  been  proposed.  On  pp.  154-155  of  Prenter  are  cited  more  than 
a dozen  papers  with  proposed  definitions. 

3.  Why  use  a finite  element  method?  The  equations,  like  (1.13), 

to  be  solved  for  the  c^'s  are  similar  to  the  equations  which  are  to  be 

solved  in  the  finite  difference  approach.  Also,  the  c^'s  that  are  obtained 

are  only  approximations  for  u at  the  grid  points.  So  the  results,  namely 

approximations  to  u at  a set  of  grid  points,  are  similar  to  those  obtained 

I 

by  the  finite  difference  approach.  Nor  does  the  finite  element  approach 
seem  to  produce  appreciably  smaller  errors  (or  appreciably  larger,  either). 
However,  there  do  seem  to  be  two  situations  in  which  the  finite  element 
method  appears  preferable. 

The  first  is  the  case  where  S2  has  curved  boundaries.  Thus,  let 
SI  be  a circle,  as  in  Figure  2.  For  the  28  triangles  that  border  on  the 
circle,  we  have  the  mismatch  between  the  straight  sides  of  the  triangles 
and  the  circle.  However,  these  triangles  come  fairly  close  to  matching 
the  circle.  It  certainly  would  not  give  a very  bad  approximation  if  we 
should  solve  (1.13)  with  the  grid  points  shown  in  Figure  2. 

A possible  improvement  would  be  to  cut  each  of  the  28  triangles 
that  border  on  the  circle  in  half  by  a line  from  the  inward  vertex  to  the 
* circle.  This  would  give  us  28  more  grid  points,  but  each  would  be  on  the 
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circle.  If  values  for  u are  assigned  on  the  boundary,  these  extra  grid 

points  would  cause  no  difficulty  at  all.  Wa  would  now  have  56  triangles 

bordering  on  the  circle,  which  would  match  it  much  better.  However, 

this  would  get  us  into  smaller  angles  for  the  triangles,  which  would  not 

2 

be  good  in  view  of  the  sin  6 in  the  denominator  on  the  right  of  (Z.  3). 

Another  way  would  be  to  add  also  new  grid  points  at  the  mid- 
points of  sides  of  triangles  that  do  not  border  the  circle.  The  triangles 


Figure  3. 

that  do  not  border  the  circle  would  be  dissected  as  shown 


in  Figure  3,  while  those  that  border  on  the  circle  would  be  replaced  by  the 
configuration  shown  in  Figure  4.  This  would  approximately  halve  the  h 
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in  (2.  3)  while  maintaining  the  0 about  the  same.  Of  course  one  does 


this  at  the  cost  of  approximately  doubling  the  number  of  equations  to  be 
solved.  However,  we  now  have  56  triangles  bordering  on  the  circle, 
which  will  certainly  give  a good  fit.  The  problem  of  trying  to  handle 
this  by  finite  differences  seems  to  involve  more  extensive  calculations. 

A very  sketchy  discussion  is  given  on  pp.  155-167  of  Prenter  of  the 
situation  with  curved  boundaries.  References  are  given  to  some  papers  in 
which  £>  is  divided  into  pieces  which  have  some  curved  boundaries;  see 
the  remarks  about  "isoparametric  transformations.  " This  enables  one  to 
get  better  fits  along  the  curved  boundaries  of  Q.  All  in  all,  if  £2  has 
curved  boundaries,  use  of  finite  element  methods  may  well  be  quite 
advantageous . 

When  one  has  a reentrant  angle  in  r (ft),  for  instance  if  £2  is 
an  L- shaped  region,  the  first  derivative  of  u is  likely  to  be  infinite  in 
the  neighborhood  of  the  reentrant  corner.  A discussion  of  this  situation  is 
given  in  "Calculation  of  Potential  in  a Sector,  Part  I,  " by  J.  Barkley  Rosser, 
MRC-Technical  Summary  Report  Number  1535,  May  197  5.  As  a result, 
the  usual  finite  difference  approximation  will  be  very  poor  indeed  near  the 
corner.  While  the  finite  difference  equations  can  be  solved  in  such  a case, 
the  resulting  values  for  grid  points  near  the  reentrant  corner  will  be  very 
poor  approximations.  If  one  uses  a finite  element  method,  with  a uniform 


-20- 


mesh,  the  approximation  for  grid  points  near  the  reentrant  corner  will  be 
poor,  and  for  similar  reasons.  However,  it  is  very  easy  to  refine  the 
mesh  near  the  reentrant  corner;  see  Figure  5.  With  the  smaller  h near 
the  corner,  one  will 


compensate  in  part  for  the  larger  derivatives.  The  local  mesh  refinement 
makes  no  trouble  at  all  in  the  finite  element  method.  It  is  possible  to  have 
a local  mesh  refinement  with  the  finite  difference  method,  but  it  involves 
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awkward  interpolations,  and  the  resulting  equations  to  be  solved  do  not 
fit  into  any  of  the  schemes,  such  as  S.O.R.  orA.D.I.,  for  abiidging  the 
calculations . 


» 


If  one  has  discontinuous  values  assigned  for  u along  the  boundary, 
there  will  be  infinite  derivatives  of  u in  the  neighborhood  of  the  dis- 
continuities; see  "Effect  of  Discontinuous  Boundary  Conditions  on  Finite- 
difference  Solutions,"  by  J.  Barkley  Rosser,  MRC -Technical  Summary 
Report  Number  138  3,  June  197  5.  As  with  the  reentrant  corner,  the  finite 
difference  approximation  will  give  poor  approximations  near  the  discon- 
tinuities. The  same  would  be  true  for  a finite  element  solution  of  uniform 
mesh.  However,  it  is  very  easy  to  refine  the  mesh  near  the  discontinuities. 

In  fact,  one  might  consider  removing  the  discontinuity  by  the  methods 
given  on  pp.  221-222  of  Milne.  However,  there  are  cases  (such  as  the 
reentrant  corner)  where  local  refinement  of  the  mesh  is  very  advantageous; 
for  such  cases,  use  of  the  finite  element  method  is  worthwhile. 


4.  A numerical  example.  Let  us  approximate  the  solution  of 


(4.1) 


82u 


2 

d u 


= 2 


3x  ay 

inside  the  unit  square  0 < x < 1,  0 < y < 1,  subject  to  the  conditions 
that  u shall  be  zero  around  the  sides  of  the  square. 
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It  is  easily  verified  that  a solution  is  given  by 


(4.2)  u x(l-x) 


00 

\ 

m = 0 


3 , 3 

7T  (2rrHl) 


sinh(2m^»iry4sinh(2m.tlMl-Jllsin(2rri4l)TTXj 

sinh(2m+l)iT 


1 rom  this,  high  accuracy  solutions  will  be  calculated  for  comparison  with 
the  approximate  solutions. 

We  superpose  on  the  square  a grid  with  the  grid  points  l/ 8 unit 
apart  vertically  and  horizontally.  We  subdivide  the  square  into  triangles 
as  in  Figure  1.  We  then  proceed  to  calculate  the  coefficients  appearing  in 
(1.13).  We  could  calculate  the  partial  derivatives  of  the  u^  by  (2.2), 
but  there  is  a quicker  way.  We  note  that  by  (2.  2)  3u^/sx  is  constant 
over  each  triangle,  and  the  same  for  9U^,/&y.  So  if  we  can  determine 
the  values  of  these  at  one  point  in  a triangle,  we  know  the  values  throughout 
the  triangle.  In  Figure  6 we  show  a typical  area  of  support  for  a grid 
point,  G,  at  the  center  of  the  figure.  Along  FG,  the  value  of  u^  rises 
from  0 to  unity  in  a distance  of  1/ 8 unit.  So  9U^/ 9x  = 8 , which  holds 


E D 


Figure  6. 


throughout  both  the  triangles  I and  IV.  Along  the  line  from  G to  C, 

-8  , which  value  holds  throughout  both  the  triangles  III  and  VI. 

Along  ED  and  AB.  we  have  au,_  / ok  0,  which  value  holds  throughout 

both  the  triangles  II  and  V.  In  a similar  manner  we  find  the  value  of 

„u  / Qy  to  be:  8 throughout  both  the  triangles  V and  VI;  -8  throughout 
k 

both  the  triangles  I and  II;  0 throughout  both  the  triangles  III  and  IV. 
Thus  we  readily  calculate  the  products  on  the  left  of  (1.13)  inside  each 
triangle.  Integration  is  accomplished  by  multiplying  by  the  area  of  the 
triangle. 

Let  us  denote  the  grid  point  at  (^,  g ) by  G, .,  its  associated 
u by  u and  the  corresponding  multiplier  by  c . Since  u is 

k ij  ’ 1 

required  to  be  0 on  the  boundary,  everything  involving  uQ  will  be  0. 

So  the  left  side  of  (1.13)  reduces  to 

4Cij  ' C(i-1) j * C(i+l)j  " Ci(i-D  Ci(j+D  ' 

On  the  right  side  of  (1.13),  the  integration  is  performed  by  taking 
the  volume  of  uR  over  the  area  of  support.  But  u^  is  a pyramid  over 
the  area  of  support,  whose  volume  is  one  third  the  base,  3/64,  times 
the  altitude,  1.  As  we  have  chosen  f = 2,  we  multiply  this  by  2.  So 
(1.13)  takes  the  form 

_ J_ 

(4.3)  4c..  - C(i_1)j  - c(i+1)j  - ®i(J_1)  " ci(j  + l)  = 32 

for  each  interior  grid  point  G... 
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It  will  be  observed  that  these  are  exactly  the  same  equations 
that  would  result  from  the  finite  difference  method. 

Although  we  have  49  equations  in  49  unknowns,  one  can  greatly 
simplify  the  calculations  by  appealing  to  symmetry.  The  solution  is 
obviously  symmetric  about  each  diagonal  of  the  square,  and  about  the 
vertical  and  horizontal  lines  through  the  center.  Consider  the  case  of 
(4.  3)  that  results  when  i = j = 1.  Since  u is  zero  around  the  boundary, 
we  have  c ^ c = 0.  Also,  by  symmetry,  we  have  c ^ = c^.  So  we 

get 


(4.4) 


4c 


11 


32  ‘ 


Proceeding  similarly,  we  get 


(4.  5) 

4°2l 

C11  " C31  ' °22  = 32 

(4.6) 

4C31  ' 

C2l  ' °4l  " C32  = Jz 

(4.7) 

4C41  ' 

ZC3l  ‘ C42  = 32 

(4.8) 

4C22  “ 

2C32  " 2C2l  = 32 

(4.9) 

4C32" 

C22  “ C42  ‘ C31  ‘ °33 

(4.10) 

4C42“ 

2C32  " °41  " °43  = 32 

(4.11) 

4C33  ‘ 

2C43  ' 2C32  = *32 

(4.12) 

4c 

43 

2C33  ' °42  ” C44  “ Tz 

(4.13) 

4c  - 

44 

4c 

43  32  • 
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We  get  (4.7)  by  using  c , c which  follows  by  symmetry.  We  get 

31  51 

(4.13)  by  using  c . _ c - c c which  follows  by  symmetry. 

4 3 4 5 34  54 

These  ten  equations  are  still  quite  sparse,  so  that  their  solution 

is  quite  easy;  it  can  be  carried  out  by  a hand  calculation.  Note  that  (4.6) 

gives  c ^ in  terms  of  three  other  c's.  Now  c.^  occurs  in  only  three 

other  equations,  so  that  one  can  easily  eliminate  c„. ^ . We  proceed  to 

eliminate  c,,,  c c__,  c._,  c , and  c(jl  by  means  of  equations  (4.4), 
li  31  lc  4Z  .5  3 44 

(4.6),  (4.8),  (4.10),  (4.11),  and  (4.13).  This  is  an  instance  of  using  a 
consistent  ordering  of  the  mesh  points;  see  Smith,  pp.  85-86. 

There  results 


(4.14) 

llc  21  - 

°41  " 

(4.15) 

9C32" 

3C21 

(4.16) 

13c41  " 

2C21 

(4.17) 

?C43- 

6C32 

3c„  — 

3 2 32 


‘ 2C41  ‘ 3C43 


- 4c  - c 
32  43 


"41  4 ' 


1 

4 

_7_ 

32 


We  use  (4.17)  to  eliminate  c , from  the  other  three  equations, 

41 

after  which  the  result  of  (4.16)  can  be  used  to  eliminate  c That 
brings  us  down  to  two  equations  in  two  unknowns,  and  we  get  finally: 
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c 


6 19 
17408 


0.035558 


1 1 

966 

C 2 1 ' 17408 

1 146 
C 3 1 = 17408 

1 20  2 
C41  " 17408 


£ 0.055492 
£ 0.065832 
£ 0.069049. 


From  these,  we  get  c^,  c^>  anc*  C42  (4.6),  and  (4.7). 

Then  we  get  c„,  and  c,,  by  (4. 9 ) and  (4 . 10) . Finally  we  get  c 

33  43  44 

by  (4.13).  Approximate  values  are  listed  in  Table  1. 

We  can  improve  our  estimates  by  using  Richardson's  deferred 

approach  to  the  limit;  see  Smith,  pp.  140-141.  This  is  not  only  valid 

for  the  finite  difference  method,  but  also  for  the  finite  element  method 

2 

because  in  this  case  the  error  also  varies  as  h ; see  Prenter,  p.  253. 
For  this,  we  take  a grid  with  1/4  unit  spacing  between  the  grid  points. 
This  gives 


_ _H_ 
C22  " 1 28 

14 

°42  128 


C44  1 28 


If  we  extrapolate  from  these,  we  get  the  values  listed  in  Table  1 . 
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TOP  ENTRY  IS  ACCURATE 

SECOND  ENTRY  FROM  (4.4)  TO  (4.13) 

THIRD  ENTRY  IS  RICHARDSON  EXTRAPOLATION 
FOURTH  ENTRY  BASED  ON  (4.  21) 


0 . 1 47  3 
0.  1456 
0 . 1 47  2 
0.1405 


0. 1321 
0. 1305 
0. 1323 
0.  1 287 


0.  1 394 
0. 1377 
0.  1394 
0. 1353 


0.0906  0. 

0.0893  0. 

0.0905  0. 

0.0886  0. 


1089 

0. 

1 147 

107  5 

0. 

1 133 

1089 

0. 

1 146 

106  5 

0. 

1 1 20 

0.0364 
0.0356 
0. 0 38  2 
0.0354 


0.0564  0.0667  0.0699 
0.0555  0.0658  0.0690 
0.0567  0.0669  0.0699 
0.0551  0.0653  0.0685 


Table  1. 
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From  the  more  accurate  extrapolated  values,  we  can  get  better 

o 

values  for  c^,  c^,  and  c^.  If  we  rotate  the  square  45  , equation 

(4 . 1)  is  invariant.  So  take  a grid  consisting  of  c(1,  c3J,  c33,  c 22,  c 42> 

and  c , and  their  reflections  about  the  lines  of  symmetry.  By  finite 
44’ 

differences,  we  will  get 

(4.!8)  4C11  " C22  = 16 

I (4-19)  4C31'C22_C42=^ 

(4.20)  4c33  - c22  - 2c4,  - c44  = I6  . 

If  we  construct  triangles  on  the  same  grid  points,  we  will  also  get  (4.18), 

(4.19),  and  (4.  20)  by  finite  elements.  The  improved  values  of  cir  C31’ 

and  c33  from  these  equations  have  been  entered  in  Table  1. 

Finally,  we  use  (4.  5),  (4.7),  (4.9),  and  (4.12)  to  get  improved 

values  for  c„,,  c„,  c„„,  and  c._;  these  have  been  entered  in  Table  1. 

21  41  32  43 

From  Table  1,  we  see  that  the  improved  values  are  never  off  by  more 
than  2 units  in  the  third  decimal,  and  mostly  off  by  no  more  than  3 units  in 
the  fourth  decimal.  Considering  that  this  resulted  from  a hand  calculation, 
the  agreement  is  very  good  indeed.  Of  course,  we  were  able  to  take 
advantage  of  unusual  symmetry.  One  cannot  count  on  doing  as  well  usually. 
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If  one  wishes  greater  accuracy,  the  obvious  approach  is  to  take 
a finer  mesh.  However,  this  considerably  increases  the  computational 
labor.  If  a very  fine  mesh  is  taken,  one  can  abridge  the  computations  by 
appealing  to  S.O.R.  or  A.D.I.  Since  the  equations  are  the  same  for  finite 
differences  and  finite  elements,  these  techniques  are  available  for  either. 

Alternatively,  one  can  use  a more  accurate  approximating  difference 
formula.  See  p.  143  of  Smith,  or  the  two  reports:  J.  Barkley  Rosser,  "Nine- 
point  Difference  Solutions  for  Poisson's  Equation,  " MRC-Technical  Summary 
Report  Number  1523,  March  197  5:  and  J.  Barkley  Rosser,  "Finite-difference 
Solution  of  Poisson's  Equation  in  Rectangles  of  Arbitrary  Shape,  " MRC- 


Technical  Summary  Report  Number  1404,  February  1974.  There  has  been 
some  work  on  developing  finite  element  approaches  so  as  to  decrease  the 
error  without  refining  the  mesh,  but  the  theory  is  not  yet  well  advanced. 


One  could  triangulate  the  square  as  shown  in  Figure  7.  This 
would  produce  the  same  set  of  equations  as  before,  except  that  (4.13) 
would  be  replaced  by 


(4.  21) 


4c  . . - 4c  = -~ 
44  4 3 48 


For  the  solution,  we  get 


46  2 

C11 

130  56 

7 20 

C 2 1 

130  56 

853 

C31 

130  56 

894 

C41 

1 30  56 

can  be 

easily  c 

results  are  listed  in  Table  1. 

These  results  are  uniformly  lower  than  those  based  on  the  other 
triangulation,  which  was  already  uniformly  low.  If  the  reader  is  dismayed 
that  two  quite  similar  triangulations  could  give  results  as  disparate  as 
these,  let  him  recall  that  these  are  only  methods  for  obtaining  APPROXIMATE 
values  of  u. 
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